
%Data = csvread('Rawdata.csv',1,0);
Data = readtable('Rawdata_1sector.csv');

data.N_i = Data.N_i;
data.N_f = Data.N_f;
data.X_i = Data.X_i*1e5;
data.X_f = Data.X_f*1e5;
data.L_i = Data.L_i*1e3;
data.L_f = Data.L_f*1e3;

pm.kappa = (data.X_i + data.X_f)./sum((data.X_i + data.X_f));
pm.t = 0.05;
pm.I = size(Data,1); %total number of sectors
pm.nu = 3;
pm.rho = 0.738;
pm.phi = 1;
pm.L= 5e6;
pm.M = 5e6;
pm.mintolp = 0.005;

Param = readmatrix('Sigmastar_1sector.csv');

pm.sigma = 0.2;
pm.se = 0.2;       
pm.theta = 1.2;
pm.bp = 1;
pm.E_i = 0.5;
pm.E_r = 50;
pm.price = Param(11);
pm.wage = Param(12);

%% SIMULATE THE THEORY MODEL
tpts   = 15;
bcgrid = linspace(1,1.5,tpts)';

Z = zeros(tpts,2);
Ninf = zeros(tpts,1);
Lc = zeros(tpts,1);
wp = zeros(tpts,1);

parfor i = 1:tpts
    [p,w,ZMat,ZstarMat, LMat,YMat,NMat,Income,VMat,TFP,Cratio] = Counterfactuals(pm,bcgrid(i),pm.bp);
    
    Z(i,:)      = ZstarMat;
    Ninf(i)   = NMat(1)./sum(NMat);
    Lc(i)     = 1 - 1./(Cratio+1);
    wp(i)     = w./p;
end
bcnorm = bcgrid;

Z = Z./Z(1,:);


 
close all
yyaxis left
scatter(bcnorm,Ninf,75,'blue','filled','LineWidth',1.5)
xlabel('b_c')
ylabel('Frac. Informal Firms')
hold on 
yyaxis right
scatter(bcnorm,Lc,75,'red','filled','^','LineWidth',1.5)
ylabel('Frac. Contract Workers')
legend('Informal firms (Ext. Margin)','Frac. cont. workers (Int. Margin)', 'Location','northwest')
saveas(gcf,[filepath.output,'/Figure1a.png'])


close all
scatter(bcnorm,Z(:,1),75,'blue','filled','LineWidth',1.5)
xlabel('b_c')
ylabel('Threshold Productivity')
hold on 
scatter(bcnorm,Z(:,2),75,'red','filled','^','LineWidth',1.5)
text(1.1,1.05,'Formal Firms', 'FontSize',14)
text(1.3,1,'Informal Firms', 'FontSize',14)
text(1.02,0.96,'Not Operational', 'FontSize',14)
legend('Informal sector','Formal sector', 'Location','northwest')
saveas(gcf,[filepath.output,'/Figure1b.png'])

fprintf('\n \n Done: Theory \n \n');
